Dynamics of driven interfaces near isotropic percolation transition 
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We consider the dynamics and kinetic roughening of interfaces embedded in uniformly random media 
near percolation treshold. In particular, we study simple discrete "forest fire" lattice models through 
Monte Carlo simulations in two and three spatial dimensions. An interface generated in the models is 
found to display complex behavior. Away from the percolation transition, the interface is self-affine 
with asymptotic dynamics consistent with the Kardar-Parisi-Zhang universality class. However, in 
the vicinity of the percolation transition, there is a different behavior at earlier times. By scaling 
arguments we show that the global scaling exponents associated with the kinetic roughening of the 
interface can be obtained from the properties of the underlying percolation cluster. Our numerical 
results are in good agreement with theory. However, we demonstrate that at the depinning transition, 
the interface as defined in the models is no longer self-affine. Finally, we compare these results to 
those obtained from a more realistic reaction-diffusion model of slow combustion. 
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I. INTRODUCTION 

Interfaces embedded in random media have received 
a considerable amount of interest recently. Such di- 
verse phenomena as pinning of flux lines in superconduc- 
tors, dynamics of flame fronts in paper and imbibition 
all contain interfaces propagating in random media with 
quenched noise [Q. For many such cases, an equation of 
motion for the d dimensional height variable t) can 
be written in the form 

1 



v\I^h{r,t) + -X\Vh{r,t)Y + F + j^{r,h), (1) 



dh{r, t) 
dt 

where F is the driving force, and the noise term rj 
represents quenched disorder and is sufficiently short- 
ranged. 

The behavior of driven interfaces near the depinning 
transition F ^ Fc at which the interface ceases to propa- 
gate and its average velocity v approaches zero has turned 
out to be nontrivial. In particular, there are two impor- 
tant universality classes that many different models of 
interface dynamics fall into at the depinning transition, 
namely the isotropic depinning (ID) or the directed per- 
colation depinning (DPD) cases [0-^. Roughly speak- 
ing, models whose microscopic dynamics is isotropic be- 
long to the ID universality class, and those with spa- 
tial anisotropy to the anisotropic universality classes, of 
which perhaps the most common one is the DPD case [|j . 
These universality classes can be distinguished by the val- 
ues of the scaling exponents associated with the interface 



near the transition, as well as the behavior of the non- 
linear term A in the equation of motion for the interface. 
For the ID case, A is kinetically generated (A ~ v) and 
vanishes at the transition, while for the DPD case this is 
no longer true. 

In this work we report the results of extensive numeri- 
cal simulations of some simple "forest fire" lattice models 
[Q where an interface propagates in a uniformly random 
background of reactants with an average concentration 
< c < 1. This is an interesting special case of a mo- 
tion of an interface through a background medium of 
quenched noise, with the additional feature that there is 
an underlying isotropic percolation transition at some fi- 
nite density c*. Below c* , the interface becomes pinned 
due to the percolation transition, and one may expect 
novel features to arise in this class of problems, which we 
call here isotropic percolation depinning. There is little 
work on the dynamics of interfaces in such isotropic lat- 
tice models, in particular near percolation HJI]. These 
type of models are also interesting from the point of view 
of recent theoretical and experimental [H,|l0| studies 
of dynamics of slow combustion in random media. 

Our results indeed reveal interesting and complex be- 
havior in the dynamics of the interface. Above the de- 
pinning transition for c > c*, the kinetic roughening of 
the interface is found to be described asymptotically by 
the Kardar-Parisi-Zhang (KPZ) [|ll[ universality class as 
described by Eq. (1) with annealed, gaussian noise. Re- 
sults consistent with the KPZ universality class were also 
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found in the simulations of Refs. of a more realis- 
tic continuum model of slow combustion. On approach- 
ing the percolation transition of the underlying lattice, A 
seems to decrease, since the nonlinear term is kinetically 
generated in the present case. We find that in this regime, 
there is a different early-time behavior. Wc show that 
in this case the global scaling exponents characterising 
the kinetic roughening of the interface can be obtained 
by utilizing results of percolation theory. In particular, 
this means that these exponents are completely deter- 
mined by the properties of the percolation cluster, and 
the continuum description of Eq. (1) must break down. 
Furthermore, we show that the interface at c* as defined 
in the models is no longer self-affine but seems to show 
multiscaling, since roughness exponents as measured nu- 
merically from different correlation functions differ . 

The results from the discrete model are compared and 
contrasted to those obtained for a continuum phase-field 
model of slow combustion introduced and studied in Refs. 
[^-||. We find that at high concentrations well above 
c*, the two models display qualitatively similar behavior. 
However, as c — » c*, the kinetic roughening of the inter- 
face is different in the two models, in that there is no 
evidence of crossover in the continuum model. We show 
through an analytic argument that this is essentially due 
to the divergence of the width of the front in the contin- 
uum model and can be understood in the framework of 
mean field theory. 

The remainder of this paper is organized as follows. 
In Section 2 the model is introduced and characterized 
in detail. In Section 3 the results of extensive Monte 
Carlo simulations in two and three spatial dimensions 
are presented. Also, a theory to explain the observed 
crossover in the dynamics is developed in this Section. A 
comparison between the discrete and continuum models 
is carried out. Finally, in Section 4 we conclude and 
discuss our results. 



and the xz plane in 3D, respectively). The dynamics of 
the model is defined by the following set of rules: During 
one Monte Carlo time step (MCS), a burning tree ignites 
all the unburned trees in a fixed, finite region around it 
and becomes a burned tree. In this work, we consider the 
nearest neighbor (NN) and next nearest neighbor (NNN) 
FF models in 2D, and the NN model in 3D. A burned 
tree will remain as such and new trees will not be gen- 
erated during simulations, in contrast to several versions 
of this basic model that display Self-Organised Critical- 
ity (SOC) ||T^. The position of the emerging interface 
/i(r, t) at column r is defined as the location of the high- 
est burning tree or the highest burned tree, if there are 
no burning trees in that column We note that this 
definition is sufficient to make the interface single valued. 

The continuum model for which we will also present 
some new results has been introduced and studied in 
Refs. 1^-^. Briefly, the model is based on a phase-field 
appraoch, utilizing a set of coupled PDE's describing the 
evolution of a thermal diffusion field r(x, y) coupled to a 
random reactants concentration field c{x,y). The inter- 
play between thermal dissipation and reaction-diffusion 
of heat generated by combustion determines the dynam- 
ics in the model. To study front propagation, the set of 
equations is discretized on a 2D lattice and solved numer- 
ically. In analogy with the FF model, the lattice sites are 
randomly filled with reactants ( "trees" ) which "burn" ac- 
cording to the kinetics defined by the PDE's. The main 
difference with respect to CA type of models is that not 
only is the dynamics more realistic, but that the effective 
range of interactions is in part determined by local com- 
bustion dynamics. Also, the interface in the continuum 
model is not sharp, but can be defined through the local 
maximum of the temperature field T{x,y). 

III. RESULTS 



II. THE "FOREST FIRE" MODELS 

We consider the following simple "Forest Fire" (FF) 
cellular automaton models IJ] on square and simple cubic 
lattices in two and three spatial dimensions, respectively. 
The status of each lattice site can be one of the following: 
(i) an empty site, (ii) a site occupied by an unburned 
tree, (iii) a site occupied by a burning tree, and (iv) a 
site occupied by a burned tree. Initially, a fraction c 
(0 < c < 1) of the sites are occupied by a tree. The 
initial distribution of trees is uniformly random with no 
spatial correlations. In the 2D case the lattice is of length 
L in the x direction with periodic boundary conditions, 
and L' in the y direction with free boundary conditions. 
In the 3D case the lattice is of length L also in the z 
direction. Unless otherwise stated, L' ^ L Jl3[ . 

The front propagation is initiated at i = by igniting 
all the trees at the bottom of the lattice (y = in 2D 



In order to quantitatively characterize the kinetic 
roughening of the interface, we have considered the fol- 
lowing quantities First, the global width w{c,t,L) 
of the interface is defined by 



w'{c,t,L) = {[h{r,t)~hir,t)]^), (2) 

where the overbar denotes a spatial average over the 
system of size L, and brackets denote configuration aver- 
aging. Correspondingly, the local width of the interface 
Wi{c,t) can be defined as 

wUc,t) = {{[h{r,t) - {h{f,t))e?)i), (3) 

where the notation ()f now denotes spatial averaging 
over all subsystems of size £ of a system of total size L. 
For growing self-affine interfaces, both the global and lo- 
cal widths satisfy the Family- Viscek scaling relation |T7| ] 
and have asymptotic behavior given by 
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^2/3, for t < L^; 



(4) 



and correspondingly for w'^{t). The quantities (3 and 
X define growth and roughness exponents, respectively, 
and X = f^z We note that in addition to using the 
width, scaling exponents can be obtained by using the 
height-height difference correlation function 



C(r, t) = {Sh{faAo) - 5h{ro + r, to + t)Y), (5) 
with dh = h — h, in the appropriate regimes [|]. 

A. Dynamics of one dimensional interfaces 

In the FF models, the emerging "fire" sweeps through 
all the sites that are connected by the nearest or next 
nearest neighbor rule throughout the system. The front 
motion can thus be sustained only for lattices whose av- 
erage concentration c is at or beyond the percolation 
treshold of a 2D square lattice, which is known to be 
c* « 0.592746 and c* « 0.407254 for the NN and NNN 
cases, respectively Thus, the pinning of the inter- 

face below c* is a direct consequence of the static perco- 
lation transition, and we call this phenomenom isotropic 
percolation depinning (IPD) here. There are two com- 
peting length scales in the problem, namely the corre- 
lation length associated with the percolation transition 
£,{c), and the lateral correlation length of the moving in- 



terface £,\\it) that grows like t^^^ 
c* 

e(c)~(c-c*)- 



In the vicinity of 



(6) 



where v defines the (static) correlation length expo- 
nent, and v = 4/3 in 2D percolation JTot . 

In the regime where pinning effects can be neglected, 
it has been demonstrated that the quenched noise in Eq. 
(1) crosses over to thermal noise In the FF models, 
this sitation is realized well above c*, where ^(c) is essen- 
tially of the order of the lattice constant and ^(c) ^ C|l (0 
readily holds. Indeed, in this regime we find that the 
interface moves with a constant velocity, and its global 
width roughens asymptotically as given by Eq. (4), with 
/3 w 1/3 and x w 1/2 in accordance with the KPZ univer- 
sality class. For example, for the NN model at c = 0.95 
we obtain /3 = 0.33(1) for L = 20000 and x = 0.50(2) for 
a system size of i = 5000. In Fig. 1 we plot w{t) vs. 
t and the effective exponent /3cff{t) — lnw{t)/ \nt that 
shows the asymptotic KPZ behavior. We would like to 
point out that there is an initial time regime where the 
width grows according to the uncorrelated random depo- 
sition model, with w{t) ^ t^^^. For the NN model, this 
regime is long-lived only for c very close to unity ||21| . 

The asymptotic KPZ behavior for c > c* is not un- 
expected, since the velocity of the interface in the FF 
models is clearly tilt-dependent, which generates the non- 
linear term proportional to A in Eq. (1). In Fig. 2(a) 



we show the behavior of A as a function of c above c* 
for the NN case. It has been calculated numerically by 
computing the average velocity v of the interface as func- 
tion of the global tilt m, and fitting a parabole to it for 
m ^ 1 ||l|,||. The interesting result is that A displays 
nonmonotonic behavior and seems to eventually decrease 
when approching the percolation transition. In fact, we 
expect that A(c) ^ as c ^ c* because the interface is 
eventually forced to propagate in the infinite percolation 
cluster which is known to be self-similar and isotropic at 
the percolation threshold . Since there is no preferred 
growth direction at c* , the tilting of the interface should 
not affect the velocity of the interface any more. 

This diminishing of A on approaching c* means that 
the nonlinear term in Eq. (1) becomes less and less im- 
portant at early times where ^(c) ^ For the con- 
tinuum description to hold, however, the local slopes of 
the interface should also remain small. We have stud- 
ied this numerically for various values of c close to c*, 
where two things can be observed for the behavior of the 
global width w{t,L). First, the range of the late-time 
KPZ scaling regime becomes smaller in time as c* is ap- 
proached from above. Second, another regime where well 
defined power-law scaling of w{t) ^ t^ can be observed 
appears at earlier times. In Fig. 3(a) we show the be- 
havior of the global width for a system of size L ~ 20000 
at c = 0.59275 (NN model). We find that starting from 
early times, there is a scaling regime where the growth 
exponent /?* = 0.88(1). Simulations of the NNN model at 
c = 0.407 give /3* « 0.88, correspondingly. In this regime 
parts of the interface become pinned by unburned regions 
on the lattice, and the interface motion consists of large 
jumps, with large local slopes appearing. This behavior 
indicates that the interfaces may not be self-affine [T^ . 

The numerically observed crossover behavior can be 
formulated theoretically by assuming that it is induced 
by the underlying percolation transition. We write the 
following scaling form for the global width w{c, t): 



«;(c,t)=e(c)/(-). 



(7) 



where Tc denotes the crossover time to the KPZ regime, 
and the scaling function f{u) has the limits 



i'^ , ff M < 1; 
'/\ if M > 1. 



(8) 



Here (3* « 0.88 and /? = 1/3. Using Eq. (5) for ^(c) 
and assuming that Tc(c) ~ (c — c*)^'^ wc find that best 
data collapse as shown in Fig. 4, is obtained for ^ = 1.3 
and A = 1.65. Taking v = 4/3, and the dynamic ex- 
ponent z* given by the exponent dmin ~ 1.13 associated 
with the scaling of the minimum path distance [p|,p3|, 
A = vz* « 1.51. Thus, our numerical results are in good 
agreement with theory. 

At c* where ^|| <C £,{c) for long times {t <C L^''), the 
interface is pinned by clusters formed by the unoccupied 
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sites, and the quenched disorder dominates. The inter- 
face follows the "edge" of the infinite percolation cluster. 
The global roughness exponent x* can be then be directly 
deduced from the geometric properties of the percolation 
transition. In particular x* = ^±/'^\\j where iy± and 
ly^ are the perpendicular and parallel correlation length 
exponents of the critical percolation cluster, respectively. 
Since the percolation cluster in the FF model is isotropic 
and v± — 1^^ = the roughness exponent = 1 m all 



dimensions. In this case, the exponent z* 



1.13, 



which leads to (3* = x*/z* = 1/dmin ~ 0.88 in excellent 
agreement with our simulations. These results indicate 
that the continuum description of Eq. (1) must break 
down at c* for the present IPD case p^ . 

We have examined the interface roughness exponent x 
numerically by studying the interface dynamics as close 
to c* as possible. We have computed the generalized qth 
order height difference correlation functions 



G9(r,i) = ([/i(r,t) -r«', forr^^ll, (9) 

by running the simulation until the interface finally 
stops (for a finite system) and approximately traces out 
the edge of the percolation cluster. For a self-afhne in- 
terface there is only one roughness exponent, and thus 
X = Xq for 8-11 Q — 2,4,6, .... Our numerical results for 
a. L = 2000 system at c = 0.5928 (NN model) give that 
X2 = 0.54(5), X4 = 0.29(3) and xe = 0.21(2). This indi- 
cates that the interface associated with the percolation 
cluster as defined in the model is not self-affine at c*. The 
reason is most likely that the overhangs in the front edge 
of the interface that follow the percolation cluster, are 
removed. However, the scaling exponents for each higher 
order correlation function that we have calculated seem 
to be very well defined, which is an indication of multi- 
scaling similar to that seen in the longitudinal structure 
functions in the study of turbulence p5[ . 

We have also numerically verified the scaling of the 
average velocity v{c) of the interface as a function of c—c* 
near the percolation threshold (see Fig. 5). It is expected 
to vanish as 



v{c) = A(c-c*)* 



(10) 



Our data for the NN model give A « 1.14 and our 
best estimate for the velocity exponent is 6* = 0.169(5) 
(the NNN model gives 9 = 0.17(5)). In order to check 
the consistency of this result, we note that there exists 
a well-known scaling relation between 9, z*, and x*, 
namely 0| 



= {z*~x*>- 



(11) 



By using the values z* — 1.13, x* — I7 and v — 4/3 we 
obtain 9 = 0.173. This is in very good agreement with 
our data. 



B. Dynamics of two dimensional interfaces 

The 3D lattice model that we have studied is a simple 
generalization of the 2D case to a simple cubic geometry. 
We only consider the NN case here. The behavior of 
the emerging surface near the percolation threshold is 
qualitatively similar to the 2D case. In particular, in the 
long time limit for c > c* the interface roughens in time 
with the growth exponent (3 — 0.24(2) as shown in Fig. 
1(b), in excellent agreement with numerical solutions of 
the d = 2 -I- 1 KPZ equation and various discrete models 
that belong to the KPZ universality class p 26 27|]. 

Closer to c* « 0.316, we see the percolation-induced 
crossover. At c = 0.316 we find numerically that the in- 
terface roughens with a growth exponent (3* — 0.72(5) 
(Fig. 3.(b)). Again, the exponents characterising the in- 
terface can be obtained from the exponents of the critical 
percolation cluster. In particular, it is reasonable to as- 
sume that the global roughness exponent x* = 1 since 
the cluster is isotropic. Moreover, the minimum path 
exponent is known in c? = 3 to be dmin = 1.38(2), and 
this determines the dynamic exponent z* = dmin 
As a consistency check, taken together with the result 
that X* = li this implies that, at depinning transition, 
/3* — X* / ~ 0.724. Our direct evaluation of j3 agrees 
very well with this prediction. We expect again that the 
interface as defined in the 3D model is not self-afhne at 
c*; however, we have not computed x* numerically. 

We have also calculated the velocity exponent and find 
9 — 0.26(2). Using the scaling exponent relation 9 = 
{z* - x*)!^ with z* = 1.38, X* = 1, and = 0.88 gives 
9 ~ 0.33 which is in reasonably good agreement with our 
data. 



C. Flame front propagation in the continuum model 

We have compared the dynamics of interface of the 
lattice model to a more realistic continuum reaction- 
diffusion model of Refs. This model is a type of 
phase-field model that couples the evolution of a ther- 
mal diffusion field to a randomly distributed concentra- 
tion field of reactants. The model couples the effects of 
thermal dissipation and diffusion to heat generated by 
combustion, via an Arrhenius-activated reaction term. 
To study front propagation, the model is discretized on 
a 2D lattice and solved numerically. In analogy with 
the FF model, the lattice sites are randomly filled with 
reactants ("trees"), with an average normalized concen- 
tration of c = c{x,y). After ignition of the bottom row 
of reactants at t = 0, the heat generated will ignite other 
occupied lattice sites around it, and the local field c(x, y) 
corresponding to the sites of the "burning" reactants will 
quickly approach zero as determined by the equations. 
A single-valued interface in the model is defined by the 
maximum of the temperature field T{x^ y) for each col- 
umn X. 
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Previously, it was shown that that the kinetic roughen- 
ing of the flame fronts generated by the continuum model 
belong to the thermal KPZ universality class In the 
limit of almost uniform background density, the KPZ de- 
scription was also derived analytically from the set of 
equations for the model The main difference with 
respect to the FF lattice model was that even very close 
to the percolation treshold of the model c* w 0.20, there 
was no evidence of percolation induced crossover. Also, 
the continuum model near c* gave results that were con- 
sistent with the mean field theory of percolation, e.g. 
V « 0.5 and d « 0.5 §. 

In Fig. 2.(b) we show the behavior of the nonlinear 
coefficient A for the continuum combustion model. Simi- 
larly to the FF lattice model, we find that A approaches 
zero for c — > c*. However, unlike the the lattice model, 
no crossover behaviour is observed as c — > c* . This is ex- 
plained as follows: From the mean-field analysis of Ref. 
0, the leading front of the thermal field decays as 

rMF(x) ^e-^/'-, (12) 

where Ijj = D /vm is the thermal diffusion length defin- 
ing the range of effective interactions in the model, and 
thus also the scale of the intrinsic thickness of the in- 
terface Wint. The constants Vm and D are the mean in- 
terface velocity and thermal diffusion constant, respec- 
tively. Using the result that w,„ ^ {c — c*)°-^, we con- 
clude that Wint ~ (c — c*)~°'^. On the other hand, in the 
MF percolation transition the correlation length scales as 
^(c) ^ {c—c*)~^^^ . These results imply that the thickness 
of the interface has the same divergence as the correla- 
tion length, within which the crossover behavior should 
be observed. Thus, everything that happens on length 
scales smaller than w-mt will be smeared out. Therefore, 
due to the increasing thickness of the interface, the sec- 
ond regime at early times is never observed. 



to indicate multiscaling. This is most likely due to the 
removal of overhangs in the way the interface is defined 
in the models. 

Comparison between the lattice models and the more 
realistic model of Refs. [^-|| was made, and qualitatively 
similar behavior was found at high concentrations. In- 
terestingly, however, the two models displayed qualita- 
tively different behavior for c c*. In particular, the 
exponents compatible with the KPZ universality were 
shown to hold for all values of c studied in Refs. [^,0. 
We demonstrate that this can be understood on a ba- 
sis on the mean-field nature of the percolation transition 
exhibited by the continuum model. 

The models studied here are particularly interesting 
from the point of view of the recent experiments on 
slow combustion of paper | |l0| , p8t . In these experiments, 
asymptotic KPZ exponents were verified for the first time 
for driven interfaces. This is in complete agreement with 
all the models here well above percolation, as well as 
the DPD universality class. Near percolation, the as- 
sumption made on the basis of the earlier experiments by 
Zhang et al. has been that DPD effects dominate [Q. 
However, the most recent experiments indicate [ p8| that 
the effective short-range exponents before KPZ asymp- 
totics may not be well defined. 
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IV. SUMMARY AND DISCUSSION 

In this work we have studied the dynamics of inter- 
faces in random media through Monte Carlo simulations 
of some discrete cellular automaton models of "forest 
fires" . We find that away from the depinning transition 
induced by the isotropic percolation transition of the un- 
derlying lattice, the kinetic roughening is asymptotically 
described by the Kardar-Parisi-Zhang ||ll|| universality 
class. In the vicinity of the IPD transition, however, the 
behavior is found to be different. At the transition, the 
global roughness exponent x* and the growth exponent 
(3* are completely determined by the geometric proper- 
ties of the percolation transition, leading to the result 
that X* = Ij /?* = 1/dmin in all dimensions. We have ver- 
ified this numerically for the exponent f3* in 2D and 3D 
cases. However, by computing the roughness exponent 
of the interface from different correlation functions, we 
find that the interface is no longer self-affine, but seems 
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Figure captions 



Fig. 1. (a) The global width w{t) vs. t in 2D for 

the NN model with c = 0.95 and L = 20000. The inset 
shows the effective growth exponent /3off(i) vs. t, where 
/3eff(t) = \nw{t)/\nt. The exact KPZ value of f3 = 1/3 
is shown by the horizontal line, (b) The global width 
w{t) vs. t for the 3D NN model, with c = 0.97 and 
L X L = 200 X 200. In both cases averages were taken 
over 100 runs. Inset shows f3e«{t), with the horizontal 
value indicating the KPZ result /3 « 0.24. 

Fig. 2. (a) A vs. c in 2D with L = 2000 for the NN 
lattice model. The data was averaged over 1000 runs, 
(b) A vs. c for the continuum model, with L = 200. 

Fig. 3. (a) w{t) vs. t in the 2D NN lattice model very 
close to the percolation transition (c = 0.59275, L = 
20000). The inset shows the effective growth exponent 
/3off(i) and the horizontal line indicates the value 0.88. 
(b) w{t) vs. t in the 3D NN lattice model very close to 
the percolation transition (c = 0.312, L x L x L = 1100 x 
1100). The inset shows the effective growth exponent 
Pefi{t) VS. t and the horizontal line indicates the value 
0.72. 

Fig. 4. Crossover scaling function f{t/Tc) of the global 
width w{c,t), as defined in Eq. (6). The unsealed data 
for different concentrations (c = 0.594, 0.60, 0.605, 0.61, 
0.615, 0.62, and 0.63, from top to bottom, and L = 1000) 
is shown in the inset. The data collapse has been ob- 
tained using ly = 1.3 and A = 1.65. See text for details. 

Fig. 5. (a) Scaling of the interface velocity ?; vs. c — c* 
for the NN 2D lattice model, with L = 2000. The straight 
line shows the best fit to the data, with 9 = 0.169 in Eq. 
(9). (b)Scaling of the interface velocity w vs. c — c* for 
the NN 3D lattice model, with L x i = 100 x 100. The 
straight line shows the best fit to the data, with 9 = 0.26. 
The error bars here are smaller than the symbol sizes. 
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